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Abstract 

Presented in this paper is a new sparse linear solver methodology mo- 
tivated by multigrid principles and based around general local transfor- 
mations that diagonalize a matrix while maintaining its sparsity. These 
transformations are approximate, but the error they introduce can be 
systematically reduced. The cost of each transformation is independent 
of matrix size but dependent on the desired accuracy and a spatial er- 
ror decay rate governed by local properties of the matrix. We test our 
method by applying a single transformation to the 2D Helmholtz equation 
at various frequencies, which illustrates the success of this approach. 

1 Introduction 

Physicists develop mathematical models from physical considerations, but the 
process of solving a model isn't always related to its physics. Intermediate 
steps of a long calculation may not have physical meaning nor grant physical 
insight. An important example of this is sparse linear system solving, which is 
used to solve discretized approximations of physical problems described by time- 
independent partial differential equations. The most general-purpose sparse lin- 
ear solvers involve either direct factorization yQ, where intermediate steps con- 
tain partially factored matrices, or iterative solvers PJ, where intermediate steps 
contain approximate solutions of increasing accuracy. For some well-understood 
problems, most notably the Poisson equation, known physical properties can be 
incorporated into a linear solver, either via multigrid methods |3] or multilevel 
preconditioning 0] , leading to algorithms that are both more physical and of op- 
timal complexity. These methods operate by approximately transforming away 
local details of the physical system, leaving successively smaller but continually 
sparse "coarsened" matrix equations that each represent the physical system on 
a different length scale. The multigrid framework [S] and coarsening procedures 
jH] have been generalized into a more algebraic formalism, but their success is 
still tied to certain spectral properties of the underlying physical system. In 
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this paper we construct a more general sparse matrix solver based on a high 
accuracy limit of algebraic multigrid that does not require as input a physical 
intuition for the problem and is not restricted by spectral properties. 

The test system we use in this paper is the 2D Helmholtz equation discretized 
on a uniform grid using finite differences and defined by the five point matrix 
stencil 

1 

1 (A -4) 1 , 
1 

where A is proportional to the frequency squared. The finite difference approxi- 
mation loses accuracy as A gets larger and breaks down entirely for A > 4 due to 
inadequate sampling of oscillations, but we are interested in the matrix problem 
and not necessarily its accuracy in reproducing the continuum problem. This is 
perhaps the simplest example of a problem for which optimal linear solvers exist 
but as yet require some analytic knowledge of the solution to construct. For 
rectangular domains, the eigenfunctions are composed of sinusoidal oscillations, 
and thus fast fourier transforms can diagonalize the matrix. Given the ex- 
act analytic inverse of the Helmholtz equation, one can hierarchically compress 
and apply it in an optimal manner using the fast multipole method | . In a 
more algebraic manner, using just the knowledge that solutions have a charac- 
teristic frequency of oscillation, it is possible to construct an optimal ray-based 
multigrid scheme Algebraic methods that don't take specific account of the 
oscillatory nature of solutions currently fail to solve the problem in an optimal 
manner. For direct factorizations, the cost has been proven to be non-optimal 
for problems on a 2D grid JU|- For the parameter range of oscillatory behavior, 
< A < 8, preconditioned based on multigrid principles fail to be optimal due 
to a loss of smoothness on coarse grids and structured direct methods fail 
due to the loss of low off-diagonal rank ^2] • 

In Section [3 we describe the form of the linear solver as a succession of 
local transformations and some of their properties and governing equations. In 
Section |31 we derive an efficient method for constructing local transformations 
and apply it to our test problem. In Section 0] we further generalize the local 
transformations by changing sparsity patterns to improve accuracy. 

2 Form of the Transformation 

The basic operation of our linear solver is to start from A, an n x n real sym- 
metric sparse matrix at some stage of factorization, and apply a real symmetric 
transformation, 

X T AX = A + E, 

that leaves us with an A that has one more diagonalized row/column than A and 
a small amount of error, E. Transformations of this form are found in direct LDL 
factorization P^, where X is the identity plus a rank one matrix and E is just 
floating point roundoff error. This form can also be related to multigrid solvers, 
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if coarsening and relaxation are combined into a single invertible transformation 
jS] and if coarsening is performed only on one single small subdomain at a time. 
The X will be a dense matrix on this small subdomain and identity outside of 
it, and E on the subdomain will be substantially larger than roundoff errors. 

The benefit of multigrid, despite the large error, is that the sparsity pattern 
of A can be more controlled and the critical filling in of A during factoriza- 
tion that prevents LDL factorization from being optimal can be avoided. The 
large error E incurred at each factorization step can be negated by including 
a multilevel refinement scheme in the linear solving procedure, but the success 
of refinement is based on details of the spectral properties of the problem [5]. 
The current prescription for reducing multigrid coarsening error is simply to 
increase fill-in of the sparsity pattern of A jS], but this relates error to matrix 
fill and reduces our ability to control the sparsity pattern. Our more general 
approach is to hold the sparsity pattern of A fixed while allowing more freedom 
in the choice of X, enough to enable \\E\\ to be arbitrarily reduced, bounded 
only by machine precision. We denote this fixed sparsity, high accuracy limit of 
algebraic coarsening as perfect algebraic coarsening and denote the X matrices 
as local sparsity-preserving transformations. 

Repeated transformations take us from our initial matrix A to a final diag- 
onal form D, 

XZ---XlA X 1 ---X n = D + O(\\E\\), 
which leads to a compact representation of the inverse of A , 

A- 1 = X 1 ■ ■■X n D- 1 X% ■■■X? + 0(\\E\\). 

If || E 1 1 can be reduced sufficiently, then this multigrid-based factorization can be 
made as accurate as a direct factorization, foregoing the need for the iterative 
steps of multigrid. If we can restrict each transformation Xi to differ from 
identity only on an n-independent sized subdomain of the problem, then each 
of these n transformations can be calculated with an n-independent cost, and 
the resulting linear solver will be of optimal 0(n) complexity. 

The restriction on each X is a "local" one, which in terms of the underlying 
grid means that a transformation that decouples a node should only act on 
neighbors of that node up to at most some m th nearest neighbor. The restricted 
transformation takes the form 

\xl 1 

/ 

where the subscript 'L' refers to a local partition and 'E' to the remaining 
external partition. In order for Eq. Q to be satisfied with a small error, we 
have to enforce the condition, A\ E Xj, = A\ E , either approximately with some 
least squares approach or exactly by finding the null space of A\ E or more 
simply by further partitioning the local region into an interior and boundary 
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and further restricting Xl to 
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Our calculation of X and A may now proceed independently of the external 
partition, with some n-independent cost dependent only on n^, rig, and nj - 



and 



'A- 



the number 



the sizes of the local, boundary, and internal partitions 
of independent nonzero elements in the symmetric All. 

We must next define an error norm to be minimized by our choice of trans- 
formation. A convenient choice of norms when dealing with variable matrices 
is the Frobenius norm, ||M||^ = \jTr{M H M). Minimizing \\E\\p directly leads 
to the expression 
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with A restricted to a given sparsity pattern and X restricted to the form in 
Eq. J5J. This error norm is problematic because it is dependent on a choice of 
normalization for X to prevent such spurious solutions as {X, A) — > and to 
prevent X from becoming singular. 

An error norm that doesn't rely on normalizing X is ||A~ T i?A'~ 1 ||i?, with 
the corresponding minimization 
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where Y = X~ x is given the same local form as X. This expression is less 
appealing because it is more nonlinear than Eq. @ in that it contains 6 th order 
variable terms rather than just quartic terms. However, it is a more direct 
minimization of the error perturbation that takes us from our approximate 
inverse to the exact inverse, 

A' 1 = XA- X X T - {XA- 1 X T ){X- T EX- 1 ){XA- 1 X T ) + 0(\\X- T EX" 1 ^). 
This error norm will be used for the remainder of this paper. 



2.1 Condition of the Transformation 

Using a local transformation to remove matrix elements is only a specific appli- 
cation of a general ability to alter the values of matrix elements while preserving 
the sparsity pattern of a matrix. We can consider a transformation to be part 
of a continuous family of transformations, X(t) and A{t), that begins at t = 
as the error free identity, X(0) = / and A(0) — A, and ends at t = 1. We evolve 
from the error free transformation by following the minimum error transforma- 
tions as we continuously turn on a non-negative constraint that enforces the 
final, restricted sparsity pattern at t = 1, 

min (f error [X(t),A(t)} + t ■ f const ra m t[X(t), A(t)]) . (5) 
x(t),A(t) V - J 
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Figure 1: Convergence of error norm with steepest descent for A = 0. 

Following this defined path of transformations, the constrained error norm in 
Eq. JSJ is non-decreasing with increasing t. In order for the final transformation 
at t = 1 to have a small error, the error must be small throughout the path and 
the Jacobian of error with respect to changes of (X, A) must have an equally 
small near-null (dX, dA) component tangent to the path. Correspondingly, we 
expect the condition number of the error minimization process to be inversely 
proportional to the minimum error attainable by the transformation. 

To illustrate the ill-conditioned nature of Eq. Q , we attempt to minimize it 
by following the negative gradient for a single transformation on our Helmholtz 
test problem at A = 0. We decouple one node on the interior of the grid without 
adding or subtracting any other terms from the sparsity pattern of A and the 
local region consists of all nodes within m hops of the decoupled node. We 
start from an initial guess of Y = I and the nonzero terms of A set to the 
corresponding values of A. At each iteration, the gradient is calculated and 
the error norm is minimized in the direction of the gradient. The error norm 
for the first 1000 iterations is plotted in Fig. ^ for several values of m. Only 
the m — 1 case converges within 1000 iterations, but the expected trend of 
decreasing error and increasing condition number with increasing m is readily 
apparent. A tractable calculation of Y and A requires a more careful treatment 
of the ill-conditioned Jacobian. 

3 Linearized Approach to Local Coarsening 

Calculating and inverting the exact Jacobian of Eq. Q is impractical due to its 
size, ill-conditioning, and large null space. The symmetric error matrix, Y T EY, 
contains \riL{nL + 1) elements to be minimized and the Y and A variables 
contain (tvj + njjni) independent unknowns. For some local partitions, such as 
m > 1 in our test problem, there are more unknowns than matrix elements to 
be minimized, but the minimization is not underdetermined due to a large null 
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space. To alleviate these difficulties we separate the minimizations of Y and 
A in an approximate way that leaves us with a well-conditioned problem in Y 
whose null space can be analytically removed and a much smaller, ill-conditioned 
problem in A. 

We approximately linearize Eq. J3J by expanding Y and A in small changes, 
Y — > (Y + dY) and A — > (A + dA), and keeping only terms to first order in dY 
and dA within the norm, 



Y T (E - dA)Y - (Pj AY) T dY - dY T {P I AY) 



Because of the restricted form of Y, dY is an nj x matrix and Pj is an 
nj x til submatrix of identity. This is not the correct way to linearize Eq. 

- there are additional linear terms proportional to \\E\\ whose neglect leads 
to a linear convergence of the minimization - but this approximation leads to 
a greatly simplified solution. The Frobenius norm is invariant with respect to 
orthogonal rotations of its operand, and we choose a particularly useful rotation 
consisting of the null space Q and the spanned space Q of {Pi AY). Due to the 
Pj, the spanned space usually contains nj vectors, but it can contain less if A 
is rank deficient. If we apply the rotation to the error norm, we can write the 
norm squared as 



Q T Y T (E - dA)YQ - {Pi AY Q) T (dY Q) - (dYQ) T (P I AYQ) 
Q T Y T (E - dA)YQ - (P I AYQ) T (dYQ)" 2 
Q T Y T (E — dA)YQ 2 



The first two terms of Eq. © can be canceled with a proper choice of dY, 



(6) 



dY 



-{PiAY)- T (y t (E - dA)Y) {I + QQ T ), 



where the inverse is a pseudo-inverse. This leaves the third term to be minimized 
by dA, 

(7) 



mm 

dA 



Q 1 Y 1 (E - dA)YQ 



which is overdetermined and ill-conditioned. 

The approximately linearized Eq. © has a significant null space, corre- 
sponding to additions to dY of the form (Pi AY)~ T S (QQ T ) for any antisym- 
metric matrix S. This null space has a size of ^n/(nj — 1) for full rank A, which 
is large enough to account for Eq. being overdetermined. 

Solving Eq. Q is the most difficult and expensive step of the error minimiza- 
tion. The cost of an unstructured QR factorization of the problem is O(n^n^). 
However, the system's matrix has some structure, it is a sum of two submatriccs 
of the Kronecker product (YQ) T ® (YQ) T . There are no existing structured QR 
factorization algorithms for this kind of matrix, but the structure allows for an 
efficient construction of the normal equations, which is a sum of two x 
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submatrices of (YQQ T Y T ) ® (YQQ T Y T ). The cost of constructing and solv- 
ing the normal equations is O(n^), which is an improvement over unstructured 
QR if <C n B . In our 2D example ~ n B making the order of complexity 
equal in both methods, but the normal equations are still faster due to a smaller 
prefactor. The disadvantage of using the normal equations is the squaring of 
the condition number, which has a noticeable effect in the ill-conditioned, small 
error limit. 

3.1 Numerical Results 

We return to the test problem at A = 0, now using the linearized solution 
approach rather than following the gradient. The same local region, A sparsity 
pattern, and initial Y and A are used as in Section 12.11 The A minimization 
is performed using the normal equations which are solved using singular value 
decomposition (SVD) for testing purposes. After each iteration the solution is 
updated, Y — » (Y + adY) and A — > (A + adA), with a chosen to minimize the 
error norm. 

The SVD of Eq. JTJl, which is performed numerically on its normal equa- 
tions, is shown in Fig. [3 for the initial Y and A. An interesting feature of 
each SVD spectra is the null space of size nj, resolved in this calculation to 
single precision, 10~ 8 , relative to the largest singular value. The null space cor- 
responds to the set of local transformations that exactly preserve sparsity and 
in this case diagonal scaling of the interior block, (V, A) — > {D^ 1 Y 1 DAD). A 
change in diagonal scaling doesn't effect the chosen error norm from Eq. J3J 
and correspondingly the error term in Eq. (JZJ) is orthogonal to the null space 
within machine precision. The contribution to dA from the null space should 
be zero, and since it can be clearly distinguished from the gap in the spectrum, 
we can simply ignore the null space component. The smallest singular value of 
the rest of the spectrum shows an exponential decay with respect to m, which 
suggests an exponential decay of the minimum error according to the argument 
in Section 12.11 The limiting effects of finite precision are clearly visible in the 
vanishing of the gap between the null and spanned space for m > 9. 

The convergence of the linearized solution approach is shown in Fig. [3J Each 
calculation takes approximately m steps to converge, which signifies the success 
of our approximate inverse Jacobian in capturing the ill-conditioned aspects 
of the problem. The error exponentially decays with m as expected from the 
spectrum of Eq. (JJJl. This spatial decay of error can be related to a spatial 
decay of A to A and Y to / by taking the error to be caused by the truncation 
of some dense exact Y to I outside a local region. The relation of the decays 
can be seen in Fig. 0] where the error norm as a function of m is plotted against 
the deviations of Y and A from / and A measured by column and plotted by 
the geometric distance on the 2D grid of the associated node from the central, 
decoupled node. Since Y and A are only defined up to a diagonal scaling of the 
interior block, the rows of Y are normalized to a 2-norm of one to make them 
unique. 

We next try the method on the more interesting < A < 8 case, though 
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Figure 2: Singular values of Eq. JJJ) for A = (calculated from the normal 
equations) . 
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Figure 3: Convergence of error norm with linearized solutions for A = 0. 
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Figure 4: Spatial decay comparison, minimum error norm for different m values 
and column norms from the m = 7 minimum error solution. 
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Figure 6: The error norm and condition number of Eq. JJJ at convergence for 
m = 4. 

only the < A < 4 range needs to be tested as the matrix for A = 4 — x can 
be mapped to A = 4 + x with a diagonal scaling. The converged error norm 
for multiple values of A and m are plotted in Fig. [3J The condition number of 
the row normalized Y is less than twelve for all calculations performed. The 
decay of A^ 1 off the diagonal is exponential in geometric distance for A < 0, 
but this qualitative change in behavior from < A < 8 doesn't cause any kinks 
in the error at A = 0. We observe that the exponential decay rate of error 
with m is approximately proportional to |A — 4|. Near the A = 4 point, the 
exponential error decay appears to break down leaving an error with an m- 
dependence proportional to the logarithmic decay of off-diagonal elements of 
A^ 1 . The most obvious matrix property to attribute to the loss of decay near 
A = 4 is the vanishing of the diagonal elements of A. 

The inverse proportionality between the condition number of the non-null 
subspace of Eq. J7J and the minimum error norm continues to hold as a function 
of A as shown in Fig. The condition number plotted is calculated at the 
converged (Y, A) value, but the condition number varies very little between 
iterations and it is within a factor of two of the condition number calculated 
from the initial (Y, A) guess. 

The loss of exponential error decay as A — > 4 signifies the disappearance of 
locally removable degrees of freedom from a model restricted in form by the re- 
striction on the sparsity pattern. As A — > 4 the wavelength of oscillations in the 
Helmholtz equation approaches four times the grid spacing, a high frequency 
limit where multigrid also fails. For the multigrid approach to continue into 
this limit, the solution must be decomposed into a sum of envelope functions 
times oscillatory solutions with wavevectors in various directions {§]. This is a 
transformation from a scalar differential equation to a vector differential equa- 
tion, which can't be represented by A — > A unless the sparsity pattern of A 
is allowed to fill in somewhat. Error decay is restored for A > 4 only because 
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the discretization of the Hclmholtz equation breaks down and the correct high 
frequency oscillations are no longer present in the matrix problem. 

4 Choice of Sparsity Pattern 

For all the tests performed in Section [3] we strictly prevented fill-in in trans- 
forming from A to A, but it is only really necessary to control fill-in enough to 
preserve the scalability of the factorization. It can be beneficial to add nonzero 
matrix elements to A because that increases the number of degrees of freedom 
in the error minimization, Eq. ijljl. and can reduce the minimum error norm. 

An important reason for filling in A is to prevent the removal of a node from 
breaking the global connectivity of a problem. The simplest case of this is a 
tridiagonal matrix, which can be associated with a problem on a ID grid. If a 
node is removed from the grid without filling in the matrix, then the grid will 
be split in half. The associated transformation would have to contain all the 
response of each half of the grid on the other and cannot in general be accurately 
made local. 

One simple way to avoid changing the connectivity of a problem is to aggre- 
gate nodes together into supernodes where all the member nodes share all the 
connections of other member nodes. Once a supernode is formed, the decoupling 
of a node in the supernode won't break any connectivity as long as one node 
remains within the supernode. The supernode concept has been used before in 
sparse Gaussian elimination for efficiency reasons .13 to allow for the use of 
dense matrix operations in inner loops, but here it serves a more fundamental 
purpose. The larger the supernodes are made, the more filled in the matrix will 
be, and the error norm will have a decreasing minimum with fixed local region 
size. If the supernodes are made large enough, then Gaussian elimination steps 
can be performed without additional matrix filling before the more expensive 
algebraic coarsening procedure in a possibly more efficient hybrid approach. For 
our example on a 2D grid, the grid of nodes can be made a grid of supernodes, 
which can be interpreted as a discretization of a vector differential equation 
where the number of vector components is the size of the supernode. 

We again return to the test problem, now with a 2D grid of supernodes 
constructed by merging p x q rectangles of neighboring nodes. A local trans- 
formation is performed to remove one node from one supernode and the local 
region is chosen to include all supernodes within m hops of the removed node. 
The converged error norm for p = 2, q = 1 is plotted in Fig. [7|and the important 
difference with Fig. |S] is the error seems to continue to decay exponentially in 
m near A — 4 rather than stagnate at A = 3.5. A comparison between three dif- 
ferent supernode sizes is plotted in Fig. |S1 For similar til all errors are roughly 
the same in the non-oscillatory regime, A < 0, and at the maximally oscillatory 
point, A = 4, while the larger supernodes' errors arc smaller in the intermediate 
oscillatory regime, < A < 4. This result suggests that supernodes are use- 
ful for increasing the rate of error decay, but if the A = 4 case can indeed be 
improved, larger supernodes are required. 
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The arguments made in Ref. |5] suggest that at least an eight wave expansion 
is required for an efficient solution of the A = 4 case, which might be properly 
captured by p — 2,q — 4 or p — 3,q — 3. However, both cost and conditioning 
are a barrier to the current approach to calculating these transformations. The 
cost of solving the normal equation version of Eq. J7J for fixed local region size 
scales as 0(p a q 3 ). The conditioning of Eq. (7J remains inversely proportional 
to the minimum error norm, but the constant of proportionality is observed to 
change substantially with p and g, causing calculations to be more ill-conditioned 
with the same minimum error norms. 

5 Conclusions 

We have studied the possibility of factoring sparse matrices by means of local 
sparsity-preserving transformations with numerical tests of a single transforma- 
tion. Intermediate stages of such a factorization require transformations to be 
performed on matrices of similar sparsity but with different values of their ma- 
trix elements, which was examined here in a simple, artificial manner by varying 
the frequency of the Hclmholtz equation. Qualitatively, we expect the interme- 
diate, "coarsened" matrices to still represent the Hclmholtz equation with the 
frequency scaled to represent a change of length scale. The absence of local 
degrees of freedom for the A = 4 case in Section l3~Tl suggests that this interpre- 
tation fails when the wavelength becomes proportional to the grid spacing. To 
continue to remove local degrees of freedom beyond this frequency, it becomes 
necessary to allow the coarsened matrices to take a more general form. 

The sparse linear solver methodology presented in this paper has demon- 
strated a behavior distinct from both direct and iterative solvers. The success 
of direct solvers is dependent on the filling in of the matrix in the intermediate 
stages of factorization, which is a graph theoretic property and is controlled by 
the order in which nodes are factored. The success of iterative solvers is depen- 
dent on the condition number of the matrix and is controlled by preconditioning 
a problem to reduce the condition number. Here the determining characteristic 
of how costly it is to solve a matrix is the decay of error of a transformation 
with respect to local region size and can be controlled by changing the local 
region or sparsity pattern. Matrix fill is no longer a problem as it is strictly 
controlled, and the error decay is a local property completely independent of 
the global spectrum and conditioning of the matrix. 

There remain technical difficulties with calculating local sparsity-preserving 
transformations that must be resolved before a practical linear solver can be 
implemented with them. The most important problem is determining whether 
or not a well-conditioned process exists for calculating local transformations. 
The ill-conditioning is associated with minimizing an error norm, and a method 
based on additional criteria might precondition the process. Another impor- 
tant problem is understanding what properties of a matrix and sparsity pattern 
determine the rate of decay of error. This is needed to determine precisely 
when sparsity patterns should be changed and how they should be changed to 
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make calculations most efficient. Once the issues associated with single transfor- 
mations are resolved, there is the additional problem of choosing the ordering 
of transformations. The ordering can determine error decay rates of succes- 
sive transformations, error propagation during factorization, and the amount 
to which the process can be parallelized. The simple answer at least for the 
purpose of parallelization is to choose as many transformations as possible on 
disjoint local regions to maximize the number of concurrent calculations of local 
transformations . 
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